#include<iostream>
#include<cmath>
#include<vector>
#include "../programming1/EquationSolver.h"
#include "LUFactorization.h"
using namespace std;

class QuadraticBSpline{
private:
	int L, R;
	vector <double> coefficient;
public:
	QuadraticBSpline(const int& x0, const int& n, const vector<double>& f, const double& f0, const double& fn) : L(x0) {
		R = L + n;
		vector <double> d(n), u(n-1), l(n-1), b(n);
		for (int i = 1; i < n-1; ++ i) {
			d[i] = 6; //对角
			l[i-1] = 1;//上面
			u[i] = 1;//下面
			b[i] = 8 * f[i];
		}
		d[0] = 5, u[0] = 1, b[0] = 8 * f[0] - 2 * f0;
		d[n-1] = 5, l[n-2] = 1, b[n-1] = 8 * f[n-1] - 2 * fn;
		vector <double> t = LU(d, u, l, b); //调用LU求解
		coefficient.resize(n+2);
		for (int i = 0; i < n; ++ i) coefficient[i+1] = t[i];
		coefficient[0] = 2 * f0 - coefficient[1];
		coefficient[n+1] = 2 * fn - coefficient[n];
	}
	double getValue(const double& _x)  {
		if (_x < L || _x > R) throw "Out of Range!";
		int i = floor(_x);
		double res = 0;
		if(i-1 >= L-1 && i-1 <= R) res += coefficient[i-1-L+1] * B(2, i-1, _x);
		if(i   >= L-1 && i   <= R) res += coefficient[i  -L+1] * B(2, i  , _x);
		if(i+1 >= L-1 && i+1 <= R) res += coefficient[i+1-L+1] * B(2, i+1, _x);
		return res;
	}
};


QuadraticBSpline QuadraticBSplineInterpolation( Function & f, const int& l, const int& n) {
	vector <double> a(n);
	for(int i = 0; i < n; ++ i) a[i] = f(l+i+0.5);
	double f0 = f(l), fn = f(l+n);
	return QuadraticBSpline(l, n, a, f0, fn);
}
